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Abstract 

The product of M complex random Gaussian matrices of size TV 
has recently been studied by Akemann, Kieburg and Wei. They show 
that, for fixed M and N, the joint probability distribution for the 
squared singular values of the product matrix forms a determinantal 
point process with correlation kernel determined by certain biorthog- 
onal polynomials that can be explicitly constructed. We find that, in 
the case M = 2, the relevant biorthogonal polynomials are actually 
special cases of multiple orthogonal polynomials associated with Mac- 
donald functions (modified Bessel functions of the second kind) which 
was first introduced by Van Assche and Yakubovich. With known re- 
sults on asymptotic zero distribution of these polynomials and general 
theory on multiple orthogonal polynomial ensembles, it is then easy to 
obtain an explicit expression for the distribution of squared singular 
values for the products of two complex random Gaussian matrices in 
the limit of large matrix dimensions. 

Keywords: singular values, products of Wishart matrices, determinan- 
tal point process, multiple orthogonal polynomial ensembles, Macdonald 
functions 

1 Introduction and statement of the results 

Products of random matrices have presently attracted the most interest due 
to their important applications in statistical physics relating to disordered 
and chaotic dynamical systems |12| . and other fields beyond physics like 
MIMO (multiple-input and multiple-output) networks in telecommunication 
[S5[ 121]) etc.. The central issue of the study is to find the distributions of 
the eigenvalues or singular values for the products in different regimes and 
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various methods have been applied to perform the spectral analysis; cf. 
recent physical literatures [H El HI [H [16] and mathematical literatures 

mm- 

This note originates from a paper by Akemann, Kieburg and Wei |3j 
concerning the correlation functions for the products of M quadratic random 
matrices with complex elements and no symmetry. More precisely, let Xj, 
j = 1, ... ,M be independent complex matrices of size N with identical, 
independent Gaussian distribution 

exp [-TriXjXj)] , (1.1) 

where the superscript * stands for conjugate transpose. The interest lies in 
the singular values of the product Pm defined by 

Pm '■= XmXm-i ...Xi. (1.2) 

Note that for M = 1, it is the well-known chiral Gaussian Unitary En- 
semble. By using change of variable and Harish-Chandra-Itzykson-Zuber 
(HCIZ) integral formula, it is shown that the joint probability density func- 
tion Vj p df(xi, . . . , xn) for the squared singular values of Pm is given by (see 
[31 Equation (2.13)]) 



,(M) 

l<i<j<N 

where 



c n i] n ^~ xi \if< N 



r Mfl 
U 0,M 



,o,j-i 



(1.3) 



M 



C™) =N\HT( 



is the normalization constant and G^ M are the so-called Meijer G-functions 
which are integrals involving Gamma functions (cf. |21[ Chapter 16]). The 
formula (|1.3p is also called "one- matrix" representation in [3]. 

Since the matrix inside the determinant in (jl.3p is labeled by indices of 
the Meijer G-fucntion, it is not convenient for the computation of /c-point 
correlation function. To overcome this difficulty, they derive an alternative 
"two-matrix" formulation in the sense of [2], which is the joint probability 
density function for the squared singular values of a single matrix X\ and the 
squared singular values of the entire product matrix Pm ■ This setting leads 
to another representation of "PjpdK^i, • • • , xn) with the help of the following 
biorthogonal polynomials. 

Let p^ M \x) and q^\y) be two sequences of monic polynomials (i.e., 
with leading coefficient one) of degree j and k respectively, and satisfy the 
biorthogonality relations 
i 

w {M) (x, y)pf I) (x) q [ M) (y) dx dy = 6^/' , j, k = 0, 1, 2, • • • , 

(1.4) 



o 
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with squared norms hj = (jl) M+1 and the weight function as 

^ M K^y) = y^G^^-J^, m>i, (i.5) 



see [3, Equations (3.2) and (3.3)]. In particular, one has 

w^Hx,y) = y- 1 e-yG^(Z\^\=y- 1 e-i-y. (1.6) 

The formulas for (x) and (y) can be derived explicitly. Indeed, we 
have 

qi M \x) = {-l) k k\L k {x) (1.7) 

is the monic Laguerre polynomials associated with weight e _t , where is 
the standard Laguerre polynomials of degree k (cf. j2TJ Chapter 18]), and 

Pk w-^irww (1 ' 8) 

see [21 Equations (3.19) and (3.21)]. 

The squared singular values of Pm then forms a determinantal point 
process of the form 

Vi P d£{xi,...,x N ) = 4t det iK^ixuXj)] , (1.9) 

1 V ! 1<1,J <iV L 

where K^ 1 ^ is the correlation kernel given by 



j=0 ">j 



with as in (|1.4|) . It is shown in [3] that the two representations (|1.3p 
and (|1.9p - (|1.10p are equivalent. 

A natural question now is to establish the limiting mean distribution of 
the squared singular values for Pm- It is observed in [3] that, after proper 
scaling, the macroscopic limits of the correlation kernel K^ 1 ^ exist for each 
M. For M = 1, it is given by the famous Marchenko-Pastur density, as 
expected. For M > 2, however, it is mentioned there the explicit formula is 
difficult to derive. Alternatively, they give the algebraic equation satisfied 
by the Stieltjes transform of the limits, following the approach in [8]. The 
relevant solution of the equation for the case M = 2 is presented as well. It 
is the aim of this note to point out that for M = 2, we still have an explicit 
formula for the limit. Our main theorem is stated as follows. 
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Theorem 1.1. Let X\ and X2 be two independent complex matrices of size 
N with Gaussian distribution as in (jl.ip . // we denote by Aj, i = 1, . . . , N 
the squared singular values of X1X2 (i.e., M = 2 in (|1.2p ), then, almost 
surely, the associated scaled empirical measure 

1 N 

■= Jj^hjN^ (1-11) 
i=l 

where 5 X stands for the Dirac point mass at x, converges weakly and in 
moments to a probability measure [i on the positive real axis with density 

*H( x) = [*r h &)> * € (°>f)> (L12) 

dx ) 0, elsewhere, 

where 

h( , 3^3 (1 + y/T=y)W - (1 - yT^y)V3 

%) = ^- ^75 > 0<y<l, (1.13) 

as N — > 00. 



This theorem also implies the statement 



lim NK^\N 2 x,N 2 x) = -r-(x), x>0, 



where K N (x,y) is defined in (jl.lOp . 

The rest of this note is mainly devoted to the proof of Theorem 11.11 The 
essential issue is the observation that the biorthogonal polynomials p k in 
(|1.8p are actually certain multiple orthogonal polynomials (cf. Nikishin and 
Sorokin [2U1 Chapter 4, §3], Ismail [TTJ Chapter 23]) that have already been 
studied by the community of orthogonal polynomials. In view of the natural 
appearance of multiple orthogonal polynomials in certain random models 
from mathematical physics over the past few years, including random matrix 
theory, non- intersecting paths, etc. (cf. [HI [19] and references therein), this 
note provides another example that bridges two different groups. We end 
this paper with some concluding remarks. 



2 Proof of the main theorem 

2.1 Multiple orthogonal polynomials associated with Mac- 
donald functions 

Multiple orthogonal polynomials are polynomials of one variable which are 
defined by orthogonality relations with respect to r different measures, where 
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r > 1 is a positive integer. They can be viewed as a generalization of or- 
thogonal polynomials, which have important applications in approximation 
theory and number theory; cf. [5| l20l 125] . 

For our purpose, let us define a function p 7 by 

Py(x) = 2x 7/2 K 7 (2 v / x), x > 0, (2.1) 

where if 7 is the Macdonald function (modified Bessel function of the second 
kind; see [211 Chapter 10]) and consider two weights 

dfM\{x) = x K p^{x) dx, dfi2(x) = x K p 1+ \(x) dx, k > — 1, 7 > 0, (2.2) 

on the positive real line. For any vector (k,m) G N 2 , the type II multiple 
orthogonal polynomials p^J^ for the system of weights {111,1x2) are such 

that pi 7 ™ is a monic polynomial of degree k + m and satisfies the following 
multiple orthogonality conditions: 







PkmfrWd^x) = 0, i = 0, (2.3) 



P [ ^ } (x)x^dfi 2 (x) = 0, j = 0,l,...,m-l. (2.4) 



By taking m = k, we set 

An explicit formula for P^ 1 '^ is given by 

k 
j=0 



where 



\jj (k + l)fc-i(« + 7 + l )k-j 

see [111 Theorem 2]. Since the measures (dpi,dp2) from (|2.2p form an AT 
system (cf. [20]), it is known that all the zeros of P^'^ are simple, lie in 
(0,+oo), and satisfy the interlacing property [10] . Furthermore, it is shown 
in [27] that Pj^'^ satisfies the following four-term recurrence relation 

xP^' K \x) = pj$(z) + b k P^' K \x) + c k P^\x) + d k P^\x) (2.7) 

with recurrence coefficients 

b k = (k + k + l)(3fc + k + 2 7 ) - (as + 1)(7 - 1), 

c k = k(k + K)(k + K + -/)(3k + 2 K + ~f), (2.8) 
d k = k(k - l)(k + k - l)(k + K)(k + k + j - l)(k + k + 7). 



5 



These polynomials constitute one of few examples of multiple orthogonal 
polynomials that are not related to the classical orthogonal polynomials. 
They are first introduced by Van Assche and Yakubovich in [27], which 
solve an open problem posed by Prunikov [22] 

We now introduce a new parameter n £ N and put 

^ ^ir^ (2-9) 

Clearly, p£*\x) is a monic polynomial of degree k for each n. For each 
, we can associate the normalized counting zero measure defined by 

p(T' K )(x)=0 

k,n V / 

A measure is called the asymptotic zero distribution of {-PfcT! } if 



lim I fMPil K) )= If (2-11) 

for every bounded continuous function / on R, i.e., it is the weak limit of the 
measures V {P^^)- Here the notation lim^ means that both k,n — > oo 
with k/n — > £ > 0. 

The explicit formula for is established by Coussement, Coussement 
and Van Assche as stated in the following proposition (see [101 Theorem 
2.7]). 

Proposition 2.1. The asymptotic zero distribution of the multiple orthog- 
onal polynomials associated with Macdonald functions (|2.9[) exists and has 
the density 



d% 0, elsewhere, 



(2.12) 



where 



hf , 3^3 (1 + yT^) 1 / 3 - (1 - VT^) 1 ^ 

^ = 1^ ^ ' < y < 1. (2-13) 

Remark 2.2. In practice, one can also scale the parameters (k, 7), say, 
putting k 1 — ^ pn and 7 1— >• qn with p, g > 0. The case when n and 7 are 
fixed corresponds to taking the limits p, q — > 0, respectively. Under this gen- 
eral setting, the normalized zero counting measure ()2.10p converges weakly 
to the first component of a vector of two measures which satisfies a vector 
equilibrium problem with two external fields and the explicit formula for the 
equilibrium vector is given in terms of solutions of an algebraic equation; 
see [2H] for details. The vector equilibrium problem that characterizes v l is 
also presented in Section 3 below. 

With the above preparations, we are ready to prove Theorem ll.il 
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2.2 Proof of Theorem ITU 

A comparison of (|1.8p and (|2.5p ~ (|2.6p immediately gives that 

Pk\x) = Pk'°\ x )i ( 2 - 14 ) 

(2) 

thus, pi are multiple orthogonal polynomials with respect to two weight 
functions (2Jf (V£), 2K 1 (2 y ^x)) for x > 0. This, together with and 
Proposition 12. II implies that, as iV — > oo, the measure 

TV - v "' ~ ''' '"'- v - x 



^:=^ V 6 x/N *=v(P§:3) (2-15) 



converges weakly to z^ 1 , which agrees with the measure /x defined in Theorem 

To show the measure \x is also the almost sure weak convergence of the 
empirical measure we note that p N is actually the average char- 

acteristic polynomials with respect to the multiple orthogonal polynomial 
ensemble (|1.9p . i.e., 



p%\x)=E 



z6C, 



" N 

-i=i 

where E refers to (jl.9p with M = 2; cf. [181 Proposition 2.2]. For a deter- 
minantal point process on R that belongs to biorthogonal ensembles [7], a 
sufficient condition that the limiting zero distribution of the average charac- 
teristic polynomials coincides with the almost sure weak convergence of the 
empirical measure of this determinantal point process is recently established 
by Hardy [J5], which in particular includes multiple orthogonal polynomial 
ensembles as a special case. By [151 Theorem 1.3 and Corollary 1.5], it 
suffices to show that 

• the recurrence coefficients for p^{n?x)/n 2N = Pfo'®'(x) are bounded 
for every N, n and certain e > such that |^ — 1| < e; 

• u ( P Nn) defined in (I2TTUD converges to rA in moments, i.e., for any 

I G N, 

lim / x l dv(P^) = [ x l d^, (2.16) 

N/n->£ J ' J 

where £ G [0, 1 +e). 

These two conditions are easily verified in our case. Indeed, from (|2.7|) and 
(EEL one has 



^k 0) (-)=4^n(-) + ^t 0) (-) 

+ c W ,„Ppi n (x) + d N , n P^l n (x), (2.17) 
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with recurrence coefficients given by 

3N(N + 1) - 1 



bN,n 



7;/ 



3£ 2 , 



7? J 



3e 4 



(2.18) 



iV 3 (iV-l) 3 , 6 

«JV,n = — > ? i 



77" 



if N/n — > £. Thus, there exists e > 0, such that 



max (l&jv 

\Z-l\<e 



n 1 1 | CN,n | j 



MjV,n|) < +00. 



(2.19) 



The equations (|2T7D and (l2"7T8j) also implies that P^'°\\) = if and only 
if A is an eigenvalue of the banded Toeplitz matrix Tjy defined by 



'-N 





i 







\ 


CN,n 


bN,n 


1 






dN,n 


CN,n 


bN,n 


1 







dN,n 


CN,n 















1 


V 









bN,n) 



(2.20) 



iVxJV 

By Gershgorin circle theorem, it follows that the spectral radius pn of Tjv 
is bounded by 

SUp \b N>n \ + SUp |c?y,n| + SUp | div,n| + 1- 
N<n N<n N<n 



This, together with (|2.19|) . implies sup^y Pn < +oo. We then have that the 
support of v(Pj!f'n) 1S uniformly bounded. Since it is already known that 

u (^ > Nn ) conver g es weakly to lA, we obtain the convergence of ^(P^'n) to 
lA in moments. 

This completes the proof of our main theorem. 



3 Concluding remarks 

In this note, we have given the explicit formula for the limiting distribution 
of squared singular values for the products of two complex random Gaussian 
matrices. Our approach is based on the finite- N density derived in [3] and 
identifying the relevant biorthogonal polynomials with the special cases of 
multiple orthogonal polynomials associated with Macdonald functions. It 
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would be interesting to see whether the general multiple orthogonal polyno- 
mials associated with Macdonald functions appear in other random matrix 
models. A possible candidate might be the products of rectangular matrices 

The next challenge is to establish the microscopic limit of the correla- 
tion kernel (jl.lOp . which will lead to universality results for the products 
of Wishart random matrices. From Theorem 11.11 it is readily seen that, 
for M = 2, the density function of the limiting mean distribution vanishes 
like a squared root at the soft edge 27/4, and blows up of order x~ 2 / 3 at 
the hard edge. This observation implies the expectation to, as also pointed 
out in [3], the universal sin and Airy kernels for the bulk and soft edge of 
the spectrum, but new universality classes are required for the description 
of the local behavior at the hard edge. Due to the appearance of multiple 
orthogonal polynomials, a possible way to tackle this problem is to perform 
Deift/Zhou steepest descent analysis [13] for the associated Riemann-Hilbert 
problem [26J. To this end, it is worthwhile to mention an alternative char- 
acterization of the measure /J, given in Theorem II. li Consider the energy 
functional defined by 



// 



log: — - — rdu 1 (x)du 1 (y) + [[ log- — - — rdu 2 {x)du 2 {y) 



x 



f [ log ■ — - — -dv\ {x)dh> 2 (y) + [ V4xdui (x) 
JJ \x-y\ J 



We want to minimize this functional over all vectors of measures {vi,v 2 ) 
such that supp(z^i) C [0,oo), J dv\ = 1 and supp(z^) C (— 00, 0], J du 2 = 
1/2. It is shown by Roman and the author in [25] that the measure /i 
is the first component of the unique minimizer. We believe such kind of 
characterization will play a important role in the nonlinear steepest descent 
analysis of the Riemann-Hilbert problem. 
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